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Let  x0,  X|,  x2.  ...  be  a  bounded  sequence  of  points  some  of  which  may  be  repeated. 
The  problem  of  Rational  Hermite  Interpolation  of  type  ( m,n )  where  m  +  n  =  N  is  to  determine 
a  rational  function  Rmn(x)  =  (/( x)/V(x)  with  deg(  t/)<m  and  deg(K)<«,  which  interpolates 
an  analytic  function  /( x)  at  the  first  N  +  1  points  of  the  sequence.  If  a  point  x ,■  is  repeated 
mi  +  1  times  then  Rmn(x)  should  interpolate  /(jc)  and  its  first  m,  derivatives  at  jc,.  Hermite 
solved  this  problem  for  (m,»)  =  (N, 0)  by  constructing  the  Hermite  Interpolating  Polynomial 
PN(x)  such  that 

N 

f(x)-PN(x)  =  (*-*,) 

1-0 

where  g(jc)  is  analytic.  The  general  problem  of  Rational  Hermite  Interpolation  is  to  find  all  R 

mn 

satisfying  m  +  n  =  N  which  also  interpolate  /(x)  i.e., 

N 

f(x)-Rmn(x)  =.  gU)J~[  (x-x,)  (1) 

1-0 

The  two  extreme  cases  for  this  problem  have  special  names  :  When  the  sequence  of  points  are 
distinct  it  is  called  Cauchy  Interpolation  and  when  all  the  points  are  the  same  it  is  called  Pad* 
Table. 

A  rational  function  Am„(x)  -  U(x)/V[x)  is  said  to  solve  the  Modified  Hermite 
Interpolation  Problem  if 

N 

l/(x)  3  /(x)P(x)  mod  (x-x,)  (2) 

1-0 
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Reproduction  in  whole  or  in  part  is  permitted  for  any  purpose  of  the 
United  States  government. 
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If  Rm„(x)  solves  equation  (1)  then  equation  (2)  is  automatically  satisfied.  However,  for  some 
choices  of  m  and  n  equation  ( 1 )  may  have  no  solutions,  and  in  that  case  there  is  a  paramet¬ 
erized  family  of  solutions  to  equation  (2).  However,  each  solution  (l/(x),F(x))  to  equation 
(2)  then  yields  the  same  rational  function.  This  unique  function  is  called  the  (m,/r)*h  Rational 
Interpolant  for  /(x).  Thus  the  set  of  rational  interpolants  for  /(x),  which  is  called  the  Rational 
Interpolation  Table  for  /(x),  contains  all  solutions  to  the  problem  of  rational  Hermite  interpola¬ 
tion. 

D.  Warner  studied  this  problem  in  his  thesis  [12].  In  [13],  he  showed  all  solutions  to 

the  Modified  Hermite  Interpolation  Problem  could  be  computed  by  Kronecker’s  Algorithm  [8]. 

We  have  independently  discovered  this  and  the  result  that  Pad6  approximants  can  be  computed 

by  Euclid’s  Algorithm  prior  to  the  paper  of  McEliece  and  Shearer  [9].  Additionally,  we  have 

shown  that  Kronecker’s  Algorithm  and  the  Extended  Euclidean  Algorithm  are  virtually  the 

same.  Our  results  go  beyond  those  of  [8,9]  to  include  new  computional  techniques  as  well  as 

theoretical  unifications. 

S  N 

Let  Un  =*  II  (x-x,)  and  U.  =  2  ajc'  be  the  Hermite  interpolation  polynomial  of  /(x)  . 
u  i-o  '  1  /- 0 

The  extended  Euclidean  algorithm  applied  to  U0  and  U ,  computes  a  sequence  of  quotients  and 
remainders  according  to  the  formula  for  division: 

Uj+i  =  together  with  iterations  for  computing  the  "comultipliers": 

IPJ+,  =  ft'j-j-Q/IVj  and  Vi+l  *  K(_,-0(K(  for  i>  1,  where  initially 
W0  =.  1,  W,  -  0,  V0  =  0,  K,  =  1. 

Now,  an  important  relation  holds  for  each  i  : 

Wi  U0  +  VV\  =  ui  • 

and  the  following  results  can  be  established  for  the  Rational  Interpolation  Table. 

Lemma  1:  Each  step  of  the  extended  Euclidean  computation  gives  rise  to  a  unique  entry  (in 
lowest  terms)  of  the  Rational  Interpolation  Table. 

I 
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Lemma  2:  The  rational  function  Ci/Vi  obtainable  via  the  extended  Euclidean  computation 
yields  deg(£>,)  equal  entries  of  the  Rational  Interpolation  Table  along  the 
(m  +  n)Ih  anti-diagonal. 

Theorem  1  (Euclid-Hermite):  All  entries  along  the  ( m  +  /»),h  anti-diagonal  of  the  Rational 
Interpolation  Table  for  the  analytic  function  /(x)  are  computed  uniquely  by  the 
extended  Euclidean  algorithm. 

Lemma  1  and  2  and  Theorem  1  have  their  Cauchy  and  Padd  counterparts.  The  Pade 

Table  is  well  known  and  has  been  extensively  studied;  see  [3]  for  an  excellent  survey  article. 

As  an  example  we  state  the  above  results  in  the  Padd  case.  Let  x0  =  xt  =  ...  =  0, 

N 

Uq(x)  *  **♦'  and  U |(x)  =  5^ ap'  be  the  first  N  +  1  terms  of  the  Maclaurin  expansion  of 
fix).  Assuming  the  usual  definition  for  the  Padd  Table,  we  have  the  following  results: 

Lemma  IP:  Each  step  of  the  extended  Euclidean  computation  gives  rise  to  a  unique  entry  (in 
lowest  terms)  of  the  Padd  Table. 

Lemma  2P:  The  rational  function  (/,/ V,  obtainable  via  the  extended  Euclidean  computation 
yields  deg (Qf)  equal  entries  of  the  Padd  Table  along  the  (m  +  n),h  anti-diagonal. 
Theorem  IP  (Euclid-Padd):  All  entries  along  the  (m  +  n),h  anti-diagonal  of  the  Padd  Table 

for  the  Maclaurin  series  of  f(x)  are  computed  uniquely  by  the  extended  Euclidean 
algorithm. 

Fast  Computation  of  an  arbitrary  iterate  of  the  Extended  Euclidean  Algorithm 

The  computational  aspects  of  the  problems  of  the  previous  section  can  be  realized  by 
an  asymptotically  fast  extended  Euclidean  algorithm.  We  have  improved  and  extended  the 
HGCD  algorithm  of  Aho,  Hopcroft,  and  Ullman  [1]  in  two  significant  ways.  First,  we  have 
developed  an  improved  HGCD  algorithm  called  EMGCD  (for  Extended  Middle  CCD). 
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EMGCD  produces  the  2  by  3  matrix  of  polynomial  entries 


/  Uj  Wj  P,  \ 


where 


v^+./  V^+.^./Vt/,/ 


The  cost  of  EMGCD  is  less  than  the  cost  of  HGCD;  however,  both  algorithms  have  an 
0(/V  log27V)  asymptotic  cost.  Thus  Uj  and  Uj+l  are  computed  free  relative  to  HGCD.  Note 
also  that  algorithm  EMGCD  computes  all  of  U,  V,  and  W  which  are  the  essential  quantities 
of  the  extended  Euclidean  algorithm.  The  second  improvement  comes  from  generalizing 
EMGCD.  We  have  developed  algorithm  PRSDC  (Polynomial  Remainder  Sequence  by  .Divide 
and  Conquer)  which  produces  any  desired  iterate  Mj  in  the  PRS  sequence  and  not  just  the 
middle  term.  The  cost  of  PRSDC  is  also  0(N  log 2 TV) 

Algorithm  PRSDC  has  many  useful  applications.  One  example  is  the  computation  of 
the  greatest  common  divisor  of  two  polynomials  A  and  B.  By  setting  t/0(x)  =  A{x)  and 
t/,(x)  =  B(x)  and  specifying  C/k+,(x)  =0  or  deg((/k)>0  we  can  compute,  using  algorithm 


PRSDC 


Uk(x)  =  GCD(A (x),B(x))  =  Wk(x)A(x)  +  Vk(x)B(x)  . 

Another  example  of  its  utility  concerns  fast  computational  algorithms  for  the  above  Theorems. 

Using  algorithm  PRSDC  we  can  compute  an  arbitrary  entry  Rmn  where  m  +  n  =  N  of  the 

N 

Rational  Interpolation  Table  starting  with  U0  =  II  (x-x()  and  (/,  =  the  Hermite  interpolation 

/=o 

polynomial  of  /(x)  through  these  N  +  l  points.  Gustavson  [4]  has  shown,  using  the  ideas  of 
Yun  [14],  that  starting  with  xt  ,  /•/>(x()  ,  j  =  0,...,m(  ,  i  =  1 . k  that  the  Hermite  Interpola¬ 

tion  polynomial  P^(x)  through  these  k  distinct  points  can  be  found  in  0(N  log 2 AO-  Combining 
these  facts  we  can  state  the  following 

Theorem  2  (Euclid-Hermite-Cauchy-Pad6):  An  arbitrary  entry  of  the  Rational  Interpolation 
Table  for  the  analytic  function  /(x)  can  be  computed  in  0{N  log  N)  where  N  is 
the  degree  of  the  relevant  Hermite  interpolating  polynomial. 
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Fast  Toeplitz  Computation 


For  the  case  m  =  n  ,  equating  coefficients  of  xn,  x"  +  l . x2n  in  the  relation  for  the 


(w,  n)  Pad6  approximant,  we  get  a  Toeplitz  system: 

<*-  <*o\  ■'Vf' 


c  -.mo 


where  the  matrix,  denoted  by  T  ,  is  Toeplitz.  The  vectors  u  =  (ti0 . un)T  and 

v  =  (vq,...,^)7^  are  the  coefficients  of  the  (n,n)  Pad6  approximant  (l/y(x),Ky(x))  .  This  fact 
and  the  above  results  suggests  that  Euclid's  algorithm  can  be  adapted  to  solve  Toeplitz  systems 
of  equations.  We  now  state  a  new  theorem  which  is  a  compaction  of  two  theorems  due  to 
Gohberg  and  Semencul  [2].  This  theorem  reveals  that  the  computation  of  v  and  un  is,  in  fact, 
crucial. 

Theorem  3  :  Let  the  Toeplitz  matrix 


be  a  bordering  of  the  Toeplitz  matrix  T  with  one  additional  row  and  column  consisting  of  all 

T  R  T 

the  same  elements  except  two.  Suppose  x  =  (x0,...^t/|+|)  and  y  =  O’„+).  -d'0)  are 
solutions  of  Tx  =  e0  and  Ty*  =  e„+,  and  suppose  x0  =  yQ  *  0  .  Then  T  is  invertible  and 
it’s  inverse  S  is  formed  according  to  the  formula 


K*0  0  •  °\  /*>>!  *  y»\  / yn+\  0  •  0  \  /*»+!  *  *l  \  i 

x\  •  •  •  U  0  •  •  *  y*  *  •  •  i(  0  •  •  • 

*  *\  0  •  °yo  x  'ynyn+\  o  •  ox(l+l/ 


(3) 


Furthermore,  suppose  x  and  yR  solve  Tx  =  e0  and  Ty*  =  eH  and  0  =  >0  0  .  Then 

T~l  »  S  is  given  by  formula  (3)  with  x„+t  and  yH+ ,  set  equal  to  zero. 
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The  formula  (3),  for  the  system  with  T  ,  was  discovered  by  Trench  III],  used  by 
Zohar  [15]  and  given  a  convolutional  setting  by  Kailath,  Viera,  and  Morf  [7J.  In  ]7)  the 
formula  (3)  for  the  system  with  T  is  shown  to  be  the  the  discrete  analog  of  the  Christoffel 
Darboux  formula.  Suppose  now  that  Det(T)^0  .  Ordinarily  we  would  solve  Tx  =  e0  to  see  if 
x0/0  .  If  x0  =  0  then  formula  (3)  is  no  longer  valid.  However,  a_t  and  a2;,+  i  can  chosen 
so  that  Det(D*0  .  Then  x0  =  Tjf  =  Det(T)/Det(T)/0  .  Thus  we  have  the  following 
stronger  result  : 

Corollary  1  :  For  solving  Tz  =  b  it  is  always  possible  to  find  x  and  y  of  formula  (3)  such  that 
*0  =  >o  *  0 

Formula  (3)  is  important  because  it  expresses  the  inverse  S  as  a  product  of  Toeplitz 
matrices.  To  solve  Tz  =  b  we  can  form  four  matrix-vector  multiplications  to  affect  z  =  Sb. 
Now  we  observe  that  the  multiplication  of  Toeplitz  matrices  and  the  vector  b  given  by 


are  precisely  the  concatenations  of  the  four  matrices  in  formula  (3)  and  clearly  correspond  to 
polynomial  multiplications.  Performing  multiplication  modulo  f"+l  via  FFT  with  appropriate 
ordering  of  the  coefficients  x(-,  y,,  and  bjt  we  can  easily  derive  the  following  result  : 

Corollary  2  :  Given  x  and  y  with  x0  =  y0  #  0,  the  cost  of  solving  Tz  =  b  by  effecting 
z  =  Sb  without  explicitly  computing  S  =  T-1  is  0(n  log  w). 

Let  U0(x)  =  xz"  and  f/|(x)  =  £o(x'  .  The  polynomial  l/,  represents  the  Toeplitz 
matrix  T  .  Now  apply  the  extended  Euclidean  algorithm  to  U0  and  f/,  .  The  following  two 
theorems  demonstrate  the  importance  of  this  computation  and  establishes  a  direct  connection 
between  the  Euclidean  algorithm  and  the  solution  of  Toeplitz  systems  of  equations. 
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Theorem  4  :  Let  {Uj,  Vj.  Wj)  be  the  iterate  of  the  extended  Euclidean  algorithm  that  com¬ 
putes  the  (n,  n)  Pad6  approximant  to  Ut  .  Then  Det(r)^0  if  and  only  if 
deg  Wj)  =  «  • 

Theorem  5  :  Let  ( Uj,  Vj,  H ')  and  (Uj+l,  Vj+\,  Wj+\)  be  ,wo  successive  extended  Euclidean 
iterates  with  deg ({/•)  =  n.  These  two  extended  Euclidean  iterates  contain  all  the 
necessary  information  to  compute  x  and  y  where  Tx  =  e0  and  Ty  =  en  . 
Furthermore,  if  x0  -  0  then  the  same  two  extended  Euclidean  iterates  contain  all 
information  needed  to  compute  Tx  =  e0  and  Ty  =  en+)  with  jc0  =  =  1  . 

The  solutions  x  and  y  can  be  expressed  as  linear  combinations  of  the  V •  and  F/+, 
polynomials.  The  term  "all  the  necessary  information"  means  that  the  constants  of  the  linear 
combinations  turn  out  to  be  natural  by-products  of  the  extended  Euclidean  algorithm.  A  partial 
explanation  of  why  Theorem  5  is  true  is  the  fact  that  the  Pad6  Table  has  many  relationships 
(Frobenius  Identities)  connecting  the  Table  entries.  The  condition  of  Theorem  5  implies  that 
the  (/i,  n)  and  (n  -  1,  n  +  1)  Pad6  approximants  are  computed  by  successive  Euclidean 
iterates.  Theorems  4  and  5  and  formula  (3)  provide  the  basis  of  another  important  application 
of  algorithm  PRSDC.  We  state  this  application  as  follows: 

Theorem  6  (Euclid-Toeplitz)  :  The  complexity  of  solving  the  Toeplitz  system  Tz  =  b  is  at  most 
0(n  log  n)  and  the  extended  Euclidean  algorithm  can  be  used  to  effect  the 
solution  with  this  complexity. 

We  have  also  established  new  complexity  results  for  banded  Toeplitz  systems.  Let  Tbc 
be  a  banded  Toeplitz  matrix  whose  semi-bandwidths  are  b  and  c  i.e.,  a0  =  ...  =  an_b_}  =  0 
and  a„+c+\  =  ...  =  a2n  =  0  .  Then  by  applying  PRSDC  to  U0(x )  =  x"+ft+’  and 
l/,(x)  =  an+cxb+c  +  ...  +  an_b  we  can  solve  Tz  =  d  in  0(n  log  n)  +  0((b  +  c)log2(*  +  c)). 
The  best  previous  result  of  0(«  log  n )  +  0((b  +  c)  )  is  due  to  Jain  16)  and  Morf  and  Kailath 
[10,  p.  269].  Theorems  4  and  5  above  are  valid  for  the  banded  case.  The  only  change  in  their 
statements  is  the  replacement  of  (/»,«)  with  (b,n)  and  («  -  l,n  +  1)  with  (6  -  l.n  +  I)  . 
Recently,  Brent  discovered  a  fast  0(n  log'n)  algorithm  to  compute  x  and  y  via  a  fast 
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continued  fraction  expansion.  A  joint  paper  by  him  and  the  authors  is  planned  to  detail  some 
of  the  results  described  here.  The  best  previous  algorithm  to  solve  Toeplitz  systems  is  the 
O (n2)  algorithm  of  Trench  [11]  corresponding  to  the  Levinson  algorithm  in  the  continuum. 

The  Berlekamp  Algorithm,  Shift  register  synthesis,  and  BCH  decoding 

Let  S(x)  =  S|jc  +  ...  +  j2 „x2"  be  a  given  syndrome  polynomial.  The  key  equation  to 
finding  the  error  location  polynomial  of  BCH  decoding  is 

(1  +  S(x))a(x)  s  «(jc)  mod  (x2n+') 

where 

o(x)  =  1  +  ^  OjX1  and  «(jc)  =  1  +  ^  upc' 

/=  t  <= l 


and  e  =  deg(o)  =  deg(w)  is  small.  Berlekatnp’s  algorithm  is  an  0(n2)  method  [5]  for  comput¬ 
ing  o(x)  and  <a(x)  .  Algorithm  PRSDC  also  solves  this  problem.  Let  U0(x)  =  x2"+ 1  and 
Ut(x)  -  1  +  S(x)  .  Then  the  iterate  (l/.,  V j,  H'-)  of  the  extended  Euclidean  algorithm  which 
computes  the  ( n,n )  Pad6  approximant  to  (/j  is  the  solution  to  the  key  equation.  Also  the 
complexity  of  this  problem  is  lowered  to  0(n  log  n). 
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